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r^ Abstract 

I 

'~~' We present a study investigating the role of mitochondrial 
fvTariability in generating noise in eukaryotic cells. Noise in 
Jseellular physiology plays an important role in many funda- 
Okiental cellular processes, including transcription, translation, 
Ostem cell differentiation and response to medication, but the 
^Specific random influences that affect these processes have yet 
to be clearly elucidated. Here we present a mechanism by 
J;_which variability in mitochondrial volume and functionality, 
^^long with cell cycle dynamics, is linked to variability in tran- 
,__^cription rate and hence has a profound effect on downstream 
^ellular processes. Our model mechanism is supported by an 
• ^appreciable volume of recent experimental evidence, and we 
^^resent the results of several new experiments with which our 
J^odel is also consistent. We find that noise due to mitochon- 
drial variability can sometimes dominate over other extrinsic 
noise sources (such as cell cycle asynchronicity) and can signifi- 
cantly affect large-scale observable properties such as cell cycle 
length and gene expression levels. We also explore two recent 
regulatory network-based models for stem cell differentiation, 
and find that extrinsic noise in transcription rate causes ap- 
preciable variability in the behaviour of these model systems. 
These results suggest that mitochondrial and transcriptional 
variability may be an important mechanism infiuencing a large 
variety of cellular processes and properties. 



Author Summary 

Cellular variability has been found to play a major role in 
diverse and important phenomena, including stem cell differ- 
entiation and drug resistance, but the sources of this variability 
have yet to be satisfactorily explained. We propose a mecha- 



nism, supported by a substantial number of recent and new ex- 
periments, by which cell-to-cell differences in both the number 
and functionality of mitochondria - the organelles responsible 
for energy production in eukaryotes - leads to variability in 
transcription rate between cells and may hence be a significant 
source of cellular noise in many downstream processes. We 
illustrate the downstream effect of mitochondrial variability 
through simulated studies of protein expression and stem cell 
differentiation, and suggest possible experimental approaches 
to further elucidate this mechanism. 



Introduction 

Stochastic influences significantly affect a multitude of pro- 
cesses in cellular biology [T}{5]. Understanding the sources of 
this randomness within and between cells is a central current 
challenge in quantitative biology. Noise has been found to af- 
fect processes including stem cell fate decisions [6] , bet-hedging 
in bacterial phenotypes [7| |H], cancer development [5], and 
responses to apoptosis-inducing factors [TOl HJ. In this pa- 
per, we consider how mitochondria may constitute a significant 
source of this cellular noise. 

Noise in cellular processes may result from sources intrinsic 
to the gene in question (those responsible for differences in the 
expression levels of genes under identical regulation in the same 
cell) or extrinsic sources (those responsible for cell-to-cell varia- 
tion in genes under identical regulation in a population) . Both 
intrinsic and extrinsic noise sources contribute to the overall 
noise observed in, for example, transcription rates and protein 
expression levels (T2j . The interplay between intrinsic and ex- 
trinsic noise can be characterised with elegant experimental 
techniques such as dual reporter measurements .3], in which 
the expression levels of two proteins are measured within cells 
and within a population, but subtleties exist in disambiguat- 
ing intrinsic and extrinsic contributions to noise levels |13) . 
Some studies have found the contribution of extrinsic factors 
to overall noise levels to be stronger in eukaryotes [HI [T5j than 
prokaryotes [3] , although others debate this interpretation [16] . 
To investigate these influences, several mathematical models 
for the emergence of intrinsic and extrinsic cellular noise have 
been introduced and explored [T^ITTH^ . In addition, recent 
studies have investigated, both experimentally and theoreti- 
cally, the architecture of extrinsic noise and its causal factors 
P^HTOl [TOl [5SH1ZJ , though substantial uncertainty surrounds 
the importance of individual contributions (such as variability 
in cell cycle stage and cellular volume) to extrinsic noise ,28]. 

Huh and Paulsson recently argued that uneven segregration 
of cellular constituents at mitosis can contribute significantly 
to cell-to-cell differences in levels of cellular components and 
proteins in a population, focusing on stochasticity in protein 
inheritance between sister cells [5U1 |3D] . We focus on a spe- 
cific instance of this phenomenon: cell-to-cell variability in the 
mitochondrial content of cells. An experimental study per- 
formed by das Neves et al. identified uneven partitioning of 
mitochondria at mitosis as being a possibly significant source 
of extrinsic noise in eukaryotes |31j . supporting recent theo- 
retical ideas [30]. Mitochondria have been found to display 
remarkably complex behaviour interwoven with cellular pro- 
cesses [52H5i] and to display significant heterogeneity within 
cells f3TJ I551 - I57] . Mitochondrial infiuences on processes in- 
cluding stem cell differentation ,38, and cell cycle progression 



have recently been observed. 

das Neves et al. [31] observe a wide spread of mitochon- 
drial masses in a population of cells, illustrating extrinsic vari- 
ability in organelle distribution. Mitochondrial functionality 
has also been observed to vary between cells [Ml [351 H^H^ ■ 
das Neves et al. also observed a link between mitochondrial 
mass and membrane potential and cellular ATP levels, and 
found transcription rate to be a function of ATP concentra- 
tion. In addition, the modulation of mitochondrial function- 
ality, through anti- and pro-oxidant treatments, was found to 
alter cell-to-cell variability in transcription rates, with anti- 
oxidants significantly reducing variability and pro-oxidants in- 
creasing variability. These results suggest that cell-to-cell het- 
erogeneity in mitochondrial mass and functionality may prop- 
agate into extrinsic noise in transcription rate, and thenceforth 
processes further downstream, but the quantitative links be- 
hind these processes remain unclear. We introduce a simple 
approach, consistent with a range of experimental observa- 
tions, that quantitatively connects all these features and pre- 
dicts the downstream physiological influence of mitochondrial 
variability. 

Shahrezaei et al. |45j have recently shown that extrinsic 
noise can influence levels of intrinsic noise, as cell-to-cell vari- 
ability in the rates of processes such as transcription and trans- 
lation affect the intrinsic dynamics of gene expression. In ad- 
dition, they provided an extension to standard stochastic sim- 
ulation techniques to allow this variability in the production 
rates of chemical species to be accurately simulated - a prob- 
lem that has been approached using different techniques in pre- 
vious studies [ini 137] . However, this theoretical study did not 
attempt to characterise the physiological causes of this extrin- 
sic noise - an important consideration in assessing the ubiquity 
and consequences of cellular noise. Our proposal that cell-to- 
cell mitochondrial variability provides a significant source of 
extrinsic noise in transcription addresses these causes, and we 
show that extrinsic noise resulting from mitochondrial vari- 
ability could significantly influence intrinsic noise in gene ex- 
pression. 

This paper will proceed as follows. We first introduce one of 
the simplest possible mathematical models for variation in mi- 
tochondrial mass and functionality during and between cell cy- 
cles, and show that it is consistent with a wide range of experi- 
mental data, both from the literature and newly reported here, 
and allows analytical treatment. Our model includes stochas- 
tic segregation of mitochondria at mitosis and functional dif- 
ferences in mitochondria between cells, and contains a simple 
dynamic description of the time evolution of cellular volume 
and mitochondrial mass through the cell cycle. To our knowl- 
edge it is the first model of its kind which links mitochondrial 
mass and function to the cell cycle and gene expression. Wc re- 
late mitochondrial properties to the production of ATP in the 
cell, which in turn affects transcription rates: hence, variabil- 
ity in mitochondrial properties causes downstream variability 
in transcription. Next, wc incorporate the behaviour produced 
by our model into a common framework for cellular noise, and 
show that extrinsic noise due to variation in [ATP] can have 
a profound effect on gene expression levels, dominating over 
intrinsic noise. We then demonstrate the cell physiological 
implications of energy variability by showing how mitochon- 
drial variability may affect stem cell differentation. Finally, 
we discuss how our model relates to recent work characteris- 
ing sources of extrinsic noise, and suggest experiments to allow 



more refined models. 



Model 

While the heterogeneity of mitochondria has been observed 
experimentally and connected to variability in processes like 
transcription [31 and stem cell differentiation ^3S], the mecha- 
nisms by which mitochondrial variability influences other cel- 
lular processes has not been elucidated clearly. Here, we de- 
scribe a simple model which formalises these links, and note 
that it reproduces recent experimental results concerning mi- 
tochondrial heterogeneity (and variability in connected cellular 
features) |31| . The simplicity of our model means that ana- 
lytic expressions can be derived for many quantities of inter- 
est, facilitating a more complete and intuitive understanding 
of the modelled biological connections. We will then use this 
model to investigate more specific questions regarding the links 
between mitochondrial variability and transcription rate and 
stem cell differentiation. 

The central concept behind our model is illustrated in Fig. 
[l] Individual cells are characterised by three key variables: the 
volume of the cell (w); the amount of mitochondrial mass in 
the cell (n); and the degree of mitochondrial functionality (/). 
This last quantity, /, represents a coarse-grained measure of 
the efficiency of mitochondria within a cell - a factor which 
may be affected, for example, by the levels of reactive oxygen 
species (ROS), mitochondrial membrane potential, variability 
in mitochondrial protein complex abundance, and genetic dif- 
ferences between mitochondria |48| . ATP concentration in the 
cell is modelled as a function of these three quantities, and 
transcription rate is modelled as a function of ATP concentra- 
tion [31j . Variability in cell volume, mitochondrial mass and 
mitochondrial functionality arises due to stochastic inheritance 
of these quantities at cell divisions. This variability causes cell- 
to-cell differences in ATP levels, and hence transcription rate, 
in a population of cells. 



Stochastic Partitioning at Mitosis 

In our model, cells grow deterministically (see Cellular Dy- 
namics) , undergoing mitosis when their volume reaches a cutoff 
V* . When this occurs, the cell divides in two, with mitochon- 
drial mass n split stochastically between daughter cells, with 
each unit of mass being assigned to each cell with equal proba- 
bility, and cell volume also segregated randomly (see Methods) . 
In our model, the partitioning of n and v at mitosis is uncor- 
related. We use this lack of correlation both for simplicity and 
due to experimental data (see Results) illustrating that cell 
cycle length correlates well with inherited mitochondrial mass 
and poorly with inherited cell volume, indirectly suggesting a 
lack of correlation between n and v. This model was chosen as 
the most straightforward representation of stochastic division 
of discrete elements, and is likely to represent a realistic sce- 
nario if there is no explicit biological control mechanism that 
modulates the distribution of inherited mitochondria. 

Our mitochondrial mass measure n physically represents to- 
tal mitochondrial volume. However, it will be of use when 
considering the segregation of mitochondria at mitosis to con- 
sider the cell as populated by a number of discrete 'virtual' 
mitochondria. We denote these entities as 'virtual' mitochon- 
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FIG. 1: An illustration of the model we employ for mito- 
chondrial variability. This illustration qualitatively shows the 
key components of our model. Cell growth progresses determinis- 
tically according to the variables that characterise a cell: volume, 
mitochondrial mass (illustrated here by copy number) and func- 
tionality (illustrated here by shading). At mitosis, stochastic par- 
titioning occurs and daughter cells inherit a random volume, mi- 
tochondrial mass and functionality level from a parent cell. This 
stochastic inheritance leads to a heterogeneous population. Cells 
with high mitochondrial density and functionality have higher ATP 
levels, are able to grow faster, and have higher transcription rates 
than cells with lower mitochondrial mass and functionality. The 
variances associated with stochastic partitioning, the dependence 
of ATP concentration on cellular properties, and the dependence of 
growth and transcription rates on ATP are all parameters of the 
model. 



dria due to the difficulty of regarding mitochondria as individ- 
uals given the processes of fission and fusion [4^ . The system 
as chosen is scaled so as to regard n as mitochondrial copy 
number, so that, if n is measured in /xm^, each 'virtual' mito- 
chondrion possesses a default volume of 1 /xm^ (see Methods) . 
These virtual mitochondria are the discrete elements that, in 
our model, are binomially partitioned at mitosis. We use the 
binomial picture both for simplicity and due to its agreement 
with recent data on mitochondrial partitioning [3T], but note 
that a range of mitochondrial partitioning regimes have been 
observed in the literature ,291 1101 IS] , and explore (in the Re- 



sults section) the effects of wider or narrower distributions as- 
sociated with mitochondrial partitioning. 

We consider the variable / to be the degree of functionality 
of a cell's mitochondria. The inclusion of such a term is ne- 
cessitated by several experimental observations, das Neves et 
al. show that a measure of mitochondrial functionality (mem- 
brane potential) is slowly-varying with time in a given cell, 
although there is a wide distribution of functionality within a 
population {31] . It was also found that sister cells have similar 
transcriptional noise levels compared to the bulk population: 
if stochasticity were to arise from mitochondrial mass parti- 
tioning alone, we would expect sister cells, post-mitosis, to ex- 
hibit the greatest possible variation, as subsequent cell growth 
may be expected to dampen such variability [30 • Another 
experimental observation is that populations of cells treated 
with antioxidants, which improve mitochondrial functionality, 
showed a significant drop in noise levels. These results suggest 
that an extra source of noise, functional variability between 
cells, may be responsible for increasing noise levels. 

In the absence of a more refined view of functionality, we 
assume that all changes in functionality occur at division and 
that / stays constant through the cell cycle. / changes in a 
stochastic but mean-reverting fashion at division, and both 
daughters receive the same / value (see Methods for more de- 
tail) . In this simple model, the variation that a cell experiences 
due to slow changes in mitochondrial functionality through the 
cell cycle is absorbed into stochastic changes at cell division. 
We choose this modelling protocol due to the absence of de- 
tailed data on the behaviour of mitochondrial functionality on 
timescales longer than a cell cycle, and suggest that parameter- 
ising this simple system to match the experimentally observed 
distribution of mitochondrial functionality will give a reason- 
able estimate of the population variability in this quantity. In 
'Other Models' in Supplementary Information, we discuss an- 
other picture in which we allow / to vary continuously through 
the cell cycle, and show that similar results emerge when this 
alternative model is used. 

In this study, we will consider the oxidative state of a cell 
as a key mediator of its functionality /. Recent experimental 
data has shown that treating cells with pro- or anti-oxidants 
strongly affects the statistics of transcription rate variability in 
a population \Wr . Within our model, the effects of such chem- 
ical treatments on the oxidative state of cells can straightfor- 
wardly be captured by varying the parameters associated with 
functional inheritance (see Methods). 



[ATP] and Transcription Rate 

We are interested in the time evolution of [ATP] as a po- 
tential stochastic influence on downstream processes. Ref. [31j 
found ATP levels in the cell to be proportional to mitochon- 
drial mass (n) and membrane potential (a factor that may 
be absorbed into our measure of 'mitochondrial function' /), 
motivating our choice of expression for ATP concentration: 



[ATP] 



jnf 



(1) 



In this expression, 7 is a constant of proportionality linking 
the quantities within our model to a biological ATP concen- 
tration, and the meaning of the variable / now becomes ap- 
parent as a scalar multiple relating mitochondrial density to 



[ATP]. Wc note that other choices for the form of [ATP], 
including ODEs, are possible, and explore some alternatives 
in 'Other Models' (Supplementary Information), das Neves et 
al. also show a link between the total transcription rate A in 
a cell (measured through bromo-uridine incorporation across 
the whole nucleus) and [ATP], a sigmoidal curve, which we 
approximate (see 'Parameterisation of A' in Supplementary In- 
formation) with 



A = si + S2 tan~i (s3[ATP] + 54) • 



(2) 



das Neves et al. record a change in the structure of this sig- 
moid curve in experiments where cellular chromatin is artifi- 
cally decondensed. In these situations, the sigmoidal response 
of A to [ATP] becomes a hyperbolic curve, with a sharp, con- 
tinuous increase of A with [ATP] at low [ATP]. This change 
may reflect the necessity of remodelling chromatin - a process 
that requires ATP - for the transcription process. Chromatin 
remodelling has been noted by several studies [TSl HH [H] to 
play an important role in mRNA synthesis noise and hence 
downstream noise in gene expression. Rather than attempting 
to model this influence explicitly, we use the experimentally- 
determined form for A ([ATP]) to capture the overall depen- 
dence of transcription rate (including chromatin effects) on 
[ATP]. 

To summarise, in our model, transcription rate depends sig- 
moidally on ATP concentration - a relationship elucidated and 
quantified in recent experiments |31j . ATP concentration in 
turn depends linearly on the mitochondrial mass and func- 
tionality level of a cell and also on the cell volume. Cells with 
many, highly functional mitochondria will have higher levels 
of ATP and hence higher transcription rates than those with 
smaller, less functional mitochondrial populations. 



Cellular Dynamics 

Our model for cell cycle dynamics consists of equations gov- 
erning the time evolution of the key quantities volume, mito- 
chondrial mass, and mitochondrial functionality. In the light 
of a recent study ^, and as cell cycle models often assume 
the exponential growth picture, we expect an exponential form 
for cell volume growth: i) oc vF{v,n, f). Here, F{v,n,f) is a 
function expressing the dependence of volume growth rate on 
other parameters. 

We suggest that ATP concentration ([ATP]) plays a key 
role in powering growth of the cell, so cells with higher ATP 
levels have higher growth rates associated with cell volume 
and mitochondrial mass. This link postulates that biosynthe- 
sis rates are generally, like transcription, a function of ATP 
concentration. We note that although ATP concentration has 
been suggested |31j as a possible mechanism linking mitochon- 
dria and transcription rate, and some evidence supports this 
link, it may be the case that a different factor provides the 
causal mechanism, and ATP concentration is correlated with 
this underlying factor. For example, ROS, which adversely af- 
fect many cellular processes, may be an alternative to ATP, or 
a combination of ATP and ROS levels may act to determine 
transcription rate. 

Numerous historical studies, both in HeLa cells |S3| and 
other tissue types [MHS^ have found that the density p of mi- 
tochondrial mass (also called mitochondrial volume density) 
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FIG. 2: The set of data used to parameterise our model. 

Experimental data shown in blue, fitted simulated data shown in 
red. A. Ratio of larger cell volume to smaller cell volume between 
sisters at birth. B. Ratio of larger mitochondrial mass to smaller 
mitochondrial mass between sisters at birth. C. Mean and standard 
deviation of the cell cycle length in a population of cells. D. Noise 
levels in transcription rate in (C)ontrol, (A)ntioxidant-treated and 
(P)ro-oxidant-treated populations, and between (S)ister cells. Two 
other experimental values, not pictured, that were used to parame- 
terise our model are a maximum cell volume of 2 500 ^m'' (for con- 
sistency with Ref. [52]) and a mean ATP concentration of 900 /iM 
(from Ref. [60]). 



within cells of a given tissue type is consistent between gener- 
ations and within populations. This consistency suggests that 
the time evolution of mitochondrial mass should be (a) coupled 
with the time evolution of volume and (b) of a form that allows 
damping of the inherent stochasticity at mitosis. In addition 
to these features, it is presumably reasonable to assume that 
mitochondrial growth is dependent on available [ATP] (due 
to the required protein synthesis). We suggest a model that 
captures these required dependencies and incorporates mean- 
reversion, given by the dynamic equations: 



i) = afpv 
n = Pfpv, 



(3) 
(4) 



where p — n/v. 

We note that this simple model does not distinguish between 
volume growth rates at different times in the cell cycle, but 
yields a smooth exponential growth in cell size throughout the 
cell cycle. We work in this picture for simplicity and generality, 
but note that a more sophisticated model would include a more 
detailed description of cell growth as another potential source 
of variability between cells. 

The model's dynamics result (see Methods) in a convergence 
in mitochondrial density with time to a value /3/a. 



Model Parameterisation 

Values for the parameters in our model were chosen (see 
Methods) to match a subset of experimental data, illustrated 
in Fig. [2] 



Results 

In this section, we first compare recent experimental data 
to the predictions of our model and demonstrate that a good 
agreement exists across a wide range of experiments. We next 
report new experimental results of relevance to the study of mi- 
tochondrial variability and show that these too largely agree 
with the predictions from our simple model. This set of suc- 
cessful comparisons suggests that our model is capable of pro- 
ducing quantitatively sound estimates of the levels of noise as- 
sociated with mitochondrial sources of variability. Motivated 
by these results, we next show how our model allows a quan- 
titative link to be formed between mitochondrial variability 
and variability in transcription rate in cells. We explore this 
link by investigating the predictions that our model makes 
concerning noise in models of gene expression levels, and in 
models of stem cell differentiation pathways. We find that the 
mitochondrial sources of variability from our model provide a 
substantial contribution to noise levels in mRNA and protein 
levels within the cell, and can influence stem cell differenti- 
ation in a manner that depends upon the symmetry of the 
regulatory interactions that drive differentiation. 



Our simple model is sufficient to approximate a large set 
of experimental data 

Here we list a set of comparisons between predictions from 
our model and experimental studies. Unless stated otherwise, 
we will use experimental data from the study of das Neves et 
al. [31], using the protocol 'NX' to refer to data in Fig. X of 
that study. 

Distributions of mitochondrial mass and cell volume. Our 
model gives a peaked distribution skewed towards low n val- 
ues for mitochondrial mass in the bulk population (Fig. pK), 
which is similar in form to the experimental distribution (N4b) . 
The distribution of cellular volumes in a bulk population (Fig. 
^p) is found to display the quadratic decay expected from a 
theoretical treatment of cells growing exponentially [19] . 

Weak correlation between the lengths of successive cell cycles 
in a population. Fig. [Sp shows the weak relationship between 
the cell cycle length of a parent and a daughter cell, which 
qualitatively matches experimental findings (from N4h). 

Mitochondrial mass at birth is a better predictor of cell cycle 
length than cell volume at birth. Figs. l3t) andlSfe illustrate the 
correlations between cell cycle length and a cell's birth values 
of n and v respectively. The correlation between birth mito- 
chondrial mass and cell cycle length was strong (i?^ = 0.69, 
compared to the experimental value of 0.78) compared to the 
correlation between birth cell volume and cell cycle length 
{R^ = 0.22, experimental value 0.22). The same correla- 
tion behaviour is observed in experiments (from N4e and N4f ) 
which are shown for comparison. 

Transcription rate noise with cell cycle stage. We mod- 
elled progression through the cell cycle stages by assigning 
stages according to the volume w of a cell. We assign cells 
with 0.5w* <v< 0.7?;* to d, 0.7i;* <v< 0.95t;* to S, and 
0.95w* < i; < u* to G2 stages, to approximate the proportion 
of total cell cycle length that HeLa cells are observed to spend 
in each stage [ST]. Transcription rate noise was found to stay 
relatively constant (around 0.4) when population subsets at 
different positions in the cell cycle were measured (see Fig. 



3p), as observed in experiments (NSl). 

Correlation between mitochondrial mass and cell volume. 
Our model predicts a strong correlation between cell volume 
V and mitochondrial mass n (Fig. Istl). This result contrasts 
with the weak correlation observed, using forward scatter in 
flow cytometry to measure volume, by das Neves et al. (N3a) 
(we confirmed these experimental results in this study - data 
not shown). However, many historic studies have found a much 
stronger connection between mitochondrial mass and cellular 
volume. The mitochondrial density p = n/v, also referred to 
as mitochondrial volume density, has been found to exhibit 
low standard deviation (between 0.01 and 0.15 of the mean) in 
many different mammalian tissue types |53H59j and amounts 
of nitDNA have been found to display similarly low variability 
[SU [S5] . These results contrast with the extremely high vari- 
ability in mitochondrial volume density observed by das Neves 
et al. (the noise level estimated from the data is around 0.32), 
but we note that flow cytometry data (while useful for pro- 
viding approximate orderings of cells by volume) may not be 
capable of providing the absolute volume measurements which 
are required to refute the low variability in p observed in many 
other studies. 

Distribution of transcription rate per unit volume. Fig. [3jr[ 
shows the distribution of transcription rate per unit nuclear 
volume (in our model, nuclear volume is taken as proportion to 
cell volume) in the bulk population. This result follows a sim- 
ilar peaked distribution to that found experimentally (Nla). 

Others. We also note some qualitative features of our model: 
an increase in transcription rate with ATP levels is observed 
(trivially due to the functional form of A), which is also ob- 
served experimentally (N3g). We also observe an increase in 
transcription rate per unit volume with total mitochondrial 
functionality {nf in our model), found experimentally (N3d). 
Fig. |4| shows illustrative time series of the dynamic variables 
involved in simulation of our model. 



New experimental results are also consistent with this 

model 



In Fig. [3] we also present new experimental results pertain- 
ing to our model. These new experiments were designed to 
characterise two additional features of cells in a population: 
a measure of the total level of mitochondrial function within 
cells and the modulation of cell cycle lengths by changing the 
oxidative state of the cell. The total level of mitochondrial 
function is experimentally measured using the intensity of sig- 
nal from CMXRos, a dye that stains mitochondria and ac- 
cumulates according to membrane-potential, integrated over a 
whole cell (see Methods). This signal reports on the integrated 
membrane potential across the entire cell, combining measures 
of mitochondrial mass and functionality. The population dis- 
tribution of this quantity is of interest in exploring the link 
between mitochondrial mass and functionality between cells. 

The modulation of cell cycle length with cellular oxidative 
state was investigated by observing the distribution of cell cy- 
cle lengths in a control population of cells and in populations 
of cells after anti-oxidant (dithiothretiol) or pro-oxidant (di- 
amide) treatments (see Methods). Our model incorporates 
oxidative status by modulating the mean level of mitochon- 
drial functionality, so mitochondria function more readily in 
an environment with low oxidative stress than one with high 
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FIG. 3: Our simple model is consistent with experimental probes of mitochondrial and cellular variability. Comparison 
between our model (red) and experimental data (blue), following discussion in the Main Text. Experimental data from das Neves 
et al. |31| . A. Distribution of mitochondrial mass n in an unsynchronised population of cells. B. Distribution of cell volume v in an 
unsynchronised population of cells. C. Comparison of the lengths of cell cycles between generations: Gen 1 is the parent cell, Gen 2 the 
daughter. Cell cycle lengths are only weakly correlated. D. Relationship between the ratio of mitochondrial masses at birth against ratio of 
cell cycle lengths for sister pairs. E. Relationship between the ratio of cellular volumes at birth and the ratio of cell cycle lengths for sister 
pairs, showing a weaker correlation than D. F. Transcription rate noise rjx in subsets of the population in Gi, S, and G2 phases (see Main 
Text). G. Mitochondrial mass n and cell volume v are strongly correlated in our model. Some experimental evidence is contradictory (see 
Main Text). H. Distribution of transcription rate per unit volume \/v. New experimental data (see Methods). I. Distribution of total 
mitochondrial functionality (n/ in our model, CMXRos readings from experiments). J. Mean and standard deviation of cell cycle lengths 
in (A)nti-oxidant-treated, (C)ontrol, and (P)ro-oxidant-treated populations. Experimental histograms, originally presented in arbitrary 
units, have been scaled to match the mean value of the simulated data. 



oxidative stress. As mitochondrial functionality is tied in our 
model, through growth rate, to cell cycle length, we would ex- 
pect cell cycle lengths to decrease upon anti-oxidant treatment 
and increase upon pro-oxidant treatment. 

Distribution of mitochondrial functionality. Fig. [31 shows 
the distribution of total mitochondrial functionality in a pop- 
ulation of cells. In our model, this distribution is just the 
distribution of the quantity n/, and in experiments, we mea- 
sure the total membrane potential within a cell (see Methods). 
The predicted and experimentally observed distributions share 
a skewed form with similar variances. 

Cell cycle lengths in different oxidative conditions. In Fig. 
[3p we show the mean and standard deviation of cell cycle 
lengths in a control population, and upon treatment with anti- 
and pro-oxidants (see Methods). In our simulations, these 
treatments are modelled by changing the value of fc, affect- 
ing the mean functionality of mitochondria (see Table IT]) . It is 
observed that treatment with anti-oxidants reduces cell cycle 
lengths, and treatment with pro-oxidants increases cell cycle 
lengths. In our model, this behaviour emerges from the de- 



pendence of the rate of volume growth on [^TP] , and the in- 
creased [ATP] levels resulting from mitochondria with higher 
functionality. 

Mitochondrial mass and membrane potential. We also ob- 
served a linear correlation between total mitochondrial mass 
(measured with MitoGreen) and total mitochondrial mem- 
brane potential (measured with CMXRos) in experiments per- 
formed with both dyes (see Methods and 'Mitochondrial Mem- 
brane Potential' in Supplementary Information). This linear 
correlation emerges from our model due to our representation 
of total mitochondrial functionality as the product of a func- 
tional measure / with mitochondrial mass n. The observed 
correlation provides qualitative support for this representation. 



Summary of comparisons between experimental results 
and model predictions 

It can be seen that several key experimental results require 
the inclusion of terms relating to mitochondrial variability for 



an explanation. In a situation without considerable mitochon- 
drial influence on cellular variability, it may be expected that 
variability in cell cycle position among a population of un- 
synchronised cells may be a dominant source of noise. Phys- 
ical distributions subject to such cell cycle noise would be 
expected to show a variance over around an approximately 
twofold range, as this is the maximum difference in size be- 
tween two unsynchronised cells. However, several results dis- 
play data that varies over a considerably wider range than a 
factor of two, indicating that a factor other than cell cycle 
variability may be responsible. Most straightforwardly. Figs. 
[3]^ and |3j demonstrate pronounced cell-to-cell variability in 
the mass and functionality of mitochondrial populations. The 
distribution of transcription rate in Fig. [St similarly shows a 
wide range of values. 

Figs. [3JI) andjSjil demonstrate the observed fact that mito- 
chondrial inheritance at birth is a better predictor of cell cycle 
length than volume inheritance: an effect that relies on the 
presence of mitochondrial variability and mitochondrial influ- 
ence on cellular growth. The variability in cell cycle length 
observed by modulating the oxidative state of the cell in Fig. 
[3]J suggests that a source of variability that is sensitive to ox- 
idative effects strongly affects cell cycle lengths. We believe 
that these results support the hypothesis that mitochondrial 
variability provides a signiflcant contribution to the variability 
in distributions of the cellular properties we consider. 

The correspondence between experimental data and the sim- 
ulated behaviour of our model suggests that, although we have 
chosen simple functional forms in our model, the resulting be- 
haviour is biologically relevant. However, we note here that 
our model was constructed from a phenomenological philoso- 
phy, with the intention of using experimental results to con- 
struct a plausible coarse-grained explanation for the influence 
of mitochondrial variability on extrinsic noise in general and 
transcription rate in particular. Our goal was to introduce a 
simplified but consistent mathematical summary of the data 
and to use this to motivate further experiments. To this end, 
we suggest a set of experiments in 'Potential Experiments for 
Refinement' (Supplementary Information) that would support 
or contribute to further development of this model. We also 
note that many potential refinements could be made to our 
model and suggest several other functional forms in 'Other 
Models' (Supplementary Information). 
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FIG. 4: Illustration of the dynamics of our model. Example 
time series of A (transcription rate), / (mitochondrial functionality), 
n (mitochondrial mass) and v (cell volume), as a cell grows and 
divides repeatedly in our model. 



which are partitioned binomially at mitosis (see Methods). 

In Fig. [5] the functional dependence of i^x on mitochon- 
drial variability (ct^ and tr?) is shown from simulations. These 
results show that transcription rate noise is made up of signif- 
icant contributions from both mitochondrial segregation and 
functionality. We also performed simulations where (Tu, the 
variability arising from uneven volume partitioning, was set 
to zero, and where cells were sampled at the same position in 
their cell cycle, removing different ages as a source of variabil- 
ity. As Fig. [5]shows, the removal of these sources of variability 
has little impact on the overall transcription rate noise level. 
These results lead us to conclude that mitochondrial sources 
of variability provide a strong contribution to cell-to-cell vari- 
ability in transcription rate. This argument is supported by an 
approximate analytic treatment of the sources of error in tran- 
scription rate within our model (see 'Estimating Noise Contri- 
butions' in Supplementary Information). 



Mitochondrial variability can dominate noise in mRNA 
and protein expression 



Noise in transcription rate depends on noise in 
mitochondrial segregation and functionality 

We are now in a position to explore the dependence of the 
level of noise in transcription rate on the stochasticity in mi- 
tochondrial mass and function, and subsequent stochasticity 
in [ATP]. To investigate the contribution of mitochondrial 
variability sources to transcription rate noise, we performed 
simulations of our model while varying cr? , the variance asso- 
ciated with the inheritance of mitochondrial functionality, and 
fj^, the variance associated with inheritance of mitochondrial 
mass, fj^ here gives the variance of the distribution by which 
mitochondrial mass is partitioned, and varying it under the 
assumption of binomial partitioning corresponds to changing 
the mitochondrial makeup of the cell: lower cr^ corresponds to 
more mitochondrial elements, each with smaller volume, while 
higher a^ corresponds to fewer, larger mitochondrial elements. 



Having constructed and parameterised a model for mito- 
chondrial variability and its effect on transcription in the cell, 
we now investigate the connection between these factors and 
downstream quantities: mRNA expression levels, and then 
(through further extension) protein expression levels. Noise in 
protein expression levels directly affects many cellular proper- 
ties, as this noise causes cell-to-cell differences in the functional 
machinery available to perform cellular processes. Here we will 
investigate the infiuence of the mitochondrial variability sug- 
gested by the parameterisation of our model from experimen- 
tal data on existing models for mRNA and protein expression. 
We connect our findings with the substantial existing body of 
literature on this topic in the Discussion section. 

The production of mRNA and protein within a cell is of- 
ten modelled using a master equation approach, addressing 
the probability of observing a given number of molecules at a 
given time. This analytical framework lends itself to the in- 
clusion of our results for time-varying transcription rate (see 
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FIG. 5: Variability in mitochondrial mass and function- 
ality can both contribute to noise in transcription rate. 

Effects of dianging variability in mitochondrial mass inheritance 
{(Jn) and functionality (cr/) on overall transcription rate noise rjx. 
This contour plot shows the value of rix for a given combination 
of an,(Jf. More stochasticity associated with inheritance of mito- 
chondrial properties leads to higher transcription rate noise, and 
stochasticity in both mass and functional inheritance plays an im- 
portant role in transcription rate noise. Contour lines on the bottom 
surface mark different values of ri\. The 'X' mark denotes the de- 
fault parameterisation of our model. Other contour lines show that 
this relationship remains essentially identical when variability due 
to cell cycle stage and volume inheritance is removed, suggesting 
that (j„ and of are the key sources of transcription rate noise. 
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FIG. 6: Mitochondrial variability contributes strongly to 
noise in mRNA levels. Analytic and modified Gillespie simula- 
tion results for time evolution of mRNA levels with and without mi- 
tochondrial and volume variability. Bars show the mean and stan- 
dard deviation of the corresponding distribution at a given time. 
Red (+) give simulated results without inherited variability. Black 
(*) give analytic results without inherited variability. Blue (x) 
give simulated results with mitochondrial and volume variability, 
displaying much greater variance in mRNA expression. Bars are 
slightly offset in the x-direction for clarity. The inset shows two 
example time series for both simulated cases. 



Methods). Numerically, several studies have proposed tech- 
niques for incorporating time- varying rates in chemical kinetic 
systems [46l|47]: we use Shahrezaei et al.'s modification [45] 
to the Gillespie simulation method [S3] to simulate our model 
system. This protocol allows us to investigate the relative im- 
portance of intrinsic contributions (resulting in differences in 
expression levels between identical genes within a single cell) 
and extrinsic contributions (resulting in differences in expres- 
sion between identical genes in different cells in a population) . 
Fig. m] shows the increase in mRNA expression (from a level 
of zero at the start of the simulation) from our analytic ap- 
proach incorporating changing transcription rate, and in simu- 
lations run using (see Methods) a parameter set from Raj et al. 
|16j . in two scenarios: one involving only intrinsic noise effects 
(no noise due to mitochondrial variability) and one involving 
extrinsic noise in transcription rate due to mitochondrial mass, 
functionality, and cell volume variability, of the magnitudes 
found through parameterising our model with experimental 
data. It can be seen that mitochondrial variability leads to 
a large increase in the total noise in mRNA expression levels: 
without extrinsic factors, the noise in mRNA expression at a 
given time (t = 8.3 hours) was rjm — 0.04, whereas rjm — 0.40 
with extrinsic factors. We note that the means for the intrinsic 
and extrinsic noise cases differ: this result is due to the non- 
linear dependence of transcription rate on ATP concentration, 
so that [ylTP] distributions with the same mean but differ- 
ent variances may yield transcription rate distributions with 
different means. 

We can also perform simulations on the more complicated 
system involving protein production (see 'mRNA & Protein 
Levels' in Supplementary Information). With values from Raj 
et al. |16| for protein degradation and translation rate (see 
Methods), this approach allows us to simulate dual reporter 



experiments, where the expression of two distinct but iden- 
tically regulated protein-encoding genes is measured. Each 
protein was translated from a different mRNA strand, so these 
simulations tracked four quantities: the expression levels of the 
two niRNAs and the two proteins. Simulations were performed 
on synchronised and asynchronous cells, and with cr„,CT/ set 
to their model values and set to zero. In these simulations, 
mRNA molecules and proteins were also distributed binomi- 
ally between daughter cells at mitosis (see Methods). 

Dual reporter simulations performed with the parameteri- 
sation chosen from Raj et al. |16| yield very low values for 
the magnitude of intrinsic noise. This low intrinsic noise was 
found to be due to the high copy number of proteins result- 
ing from the parameterisation. To explore noise in systems 
with lower expression levels, we lowered the copy number of 
proteins by increasing the rates of mRNA and protein degra- 
dation (see Methods). Fig. [t] shows the resulting expres- 
sion levels in two proteins with and without various sources 
of extrinsic noise, at the two different degradation rate pro- 
tocols. These results show that, in our model, mitochondrial 
variability dominates the noise in protein expression levels. 
The spread of protein levels with mitochondrial and volume 
variability is much greater than the two-fold range achieved 
through cell cycle variability alone. Fig. [7] also illustrates that 
cells with higher mitochondrial mass and functionality gener- 
ally have higher protein expression levels, though inheritance 
noise makes this correlation weaker. 

In our model, we find that energy variability arising through 
mitochondrial stochasticity is the dominant source of variabil- 
ity in transcription rate, mRNA and protein expression levels. 
However, we note that the causal factors of stochasticity in 
mRNA and protein levels within the cell are significantly more 
complicated than the simple transcriptional model presented 
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FIG. 7: Effects of mitochondrial variability dominate pro- 
tein expression variability in our model. Dual reporter simu- 
lation with different sources of noise in our protein expression sim- 
ulations. All plots except (E) are normalised so that the highest 
protein expression level in the cell population is 1. Red (diamonds) 
show results from Raj et al.'s default parameterisation [16, used 
to model transcription, translation and degradation (see Methods). 
Blue (triangles) show results from this parameter set with degrada- 
tion rates increased 100-fold. Protein levels are shown from popu- 
lation of (A) unsynchronised cells with mitochondrial and volume 
variability, (B) synchronised cells with mitochondrial and volume 
variability, (C) unsynchronised cells with no mitochondrial or vol- 
ume variability, and (D) synchronised cells with no mitochondrial or 
volume variability. (E) Mean protein expression levels in the default 
parameterisation of Raj et al. with the product of mitochondrial 
mass and function nf, in the system corresponding to (A). (F) The 
equivalent plot of (A) with translation rates independent of [ATP] . 



above. The rates of many of the processes involved in more ex- 
tended models are functions of many factors which our model 
does not include. The inclusion of these complicating terms 
rapidly makes an analytic description of the model impossi- 
ble. However, we note that stochastic simulation techniques 
may be used to explore the behaviour of complex model given 
estimates for the functional dependence of process rates on 
extrinsic variables 1451. 



We also note that several studies have observed a decrease 
in intrinsic noise at higher levels of protein expression \15\ I27j . 
We do observe such a decrease, though in the default param- 
eterisation the magnitude of this effect is very small owing to 
the consistently low intrinsic noise levels. 
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FIG. 8: Transcription rate affects the stability of model 
stem cell systems. In both diagrams, curves delineate the bound- 
ary of the attractor basin corresponding to the undifferentiated cell 
state. Red (solid) to black (dotted) lines show the basin structure 
as transcription rate A increases through the given values. (A) The 
structure of the undifferentiated attractor basin in the Huang model 
given different transcription parameters, showing the widening of 
the stable undifferentiated region at high transcription rate. (B) 
The structure of the undifferentiated attractor basin in the Chickar- 
mane model, showing a decrease in undifferentiated basin size as 
transcription rate increases. The activation-repression structure of 
both models is illustrated - in (B), external terms representing the 
activation of GATAl and X exist but are set to zero in our analysis 
to allow PU.l to be expressed under some conditions. 



Mitochondrial noise, by modulating transcription rate, 
can affect stem cell differentiation 



As an illustrative application of our model, demonstrating 
its physiological relevance, we consider how, through the ex- 
trinsic effects of [ATP] on protein levels, a link between mito- 
chondrial content and stem cell differentiation behaviour may 
arise. Differentiation dynamics in stem cells have often been 
modelled as the result of expression asymmetries in lineage 
regulation genes that interact in a regulatory network [MM7] , 
but the initial sources of this expression variability have not 
been clearly elucidated and are a topic of active debate. Here 
we show that transcription rate variability resulting from mi- 
tochondrial variability can affect the dynamics of expression 
of such control genes. Experimentally, a link between stem 
cell differentiation and mitochondria was suggested by a re- 
cent study in mouse embryonic stem cells |38], showing that 
pluripotent cells with low mitochondrial membrane potential 
had higher in vitro differentiation propensity, whereas those 
with higher membrane potential remained undifferentiated and 
formed large teratomas. 

We explore two recent models for the cell fate decision be- 
tween erythroid and myeloid cell fates directed by the cross- 
antagonistic master lineage regulators GATAl and PU.l. One 
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model, by Huang et al. [35], consists of a symmetric cou- 
pled ODE system for the expression levels of these two genes, 
including cross-repression and self-activation term (see Meth- 
ods). Another model, by Chickarmane et al. [69], contains a 
similar but asymmetric ODE model, expanded to include in- 
teractions with a postulated third species which is promoted 
by GATAl and represses PU.l. The Chickarmane et al. model 
also includes external signalling terms which may act to pro- 
mote GATAl and PU.l, and repress the third species. In 
these models, cell states are defined by the relative levels of 
expression of these genes, such that undifferentiated cells have 
comparable levels of each transcription factor, while the two 
differentiated cell types correspond to a state with high lev- 
els of one factor and low levels of the other. The interactions 
between the genes are parameterised by variables such as self- 
activation and cross-repression rates (see Methods). The phase 
space for both these models comprises three attractor basins, 
corresponding to the progenitor cell type and two differenti- 
ated cell types. 

Within the Huang model, at low protein expression levels, 
smaller perturbations are required to shift attractor basins 
than at high expression levels ~ a feature consistent across 
a large range of parameterisations. Varying the parameter- 
isation of the model (modelling differentiation-inducing sig- 
nalling) changes the structure of these basins, so that the cen- 
tral undifferentiated basin becomes more or less stable to sub- 
sequent perturbation. We vary the default parameterisation 
of the model in an attempt to assess the effect of changes in 
transcription and translation rates (see Methods). We find 
that when the parameters related to the rate of production 
of proteins are low, the central, undifferentiated state is less 
stable than when they are high (see Fig. 7A), with a smaller 
volume of phase space leading to the undifferentiated basin. 

Within the Chickarmane model, a different effect is ob- 
served. As before, we investigated the volume of phase space 
corresponding to the basin representing the undifferentiated 
state. We used a nonzero value for the external signalling 
term promoting PU.l and explored the system at different 
transcription rates (see Methods). We found that increasing 
the transcription rate led to a decrease in the range of values 
of the external interaction which supported a stable undiffer- 
entiated state (see Fig. 7B). This decrease in the stability 
of the undifferentiated state arose from a smaller volume of 
phase space leading to the undifferentiated basin as transcrip- 
tion rate increased, with more phase space occupied by the 
GATAl basin. This result contrasts with the increased sta- 
bility of the Huang model at high transcription rate, due to 
the importance of the third species (the expression of which is 
dependent on transcription rate): at high transcription rate, 
the increased strength of the combined effect of self-activating 
GATAl and production of the third species shifts the basin 
structure strongly towards GATAl. 

These results suggest that cell-to-cell variability in mito- 
chondrial mass and function may, through induced variability 
in transcription rate, have a significant effect on the stability 
of bipotent cells. If differentiation dynamics are asymmetric 
and involve an intermediate species (as in the Chickarmane 
model), we find that high transcription rate destabilises the 
undifferentiated state. This destabilisation may be viewed as 
a result of the increased sensitivity of the system to pertur- 
bations: the asymmetric regulatory architecture means that a 
small increase in GATAl will be quickly amplified at high tran- 



scription rate, as more GATAl and X are quickly produced. If 
differentiation dynamics are symmetric and do not involve an- 
other species (as in the Huang model), high transcription rates 
increase the width of the basin corresponding to the undiffer- 
entiated state, acting to stabilise this state. This stabilisation 
is due to the increased robustness to perturbations afforded by 
the high production rate of both species at high transcription 
rate: without asymmetric interactions, the higher expression 
level of both genes makes the system less responsive to small 
perturbations. The results that emerge from this symmetric 
case gives results that are qualitatively comparable to an ex- 
perimental study |38| in which more cells with higher total 
mitochondrial membrane potential remained undifferentiated, 
suggesting that high mitochondrial performance stabilises the 
undifferentiated state. 

Another, higher-order effect may conceivably play a role in 
both situations: several studies have found that, at high pro- 
tein abundance levels (which may result from high transcrip- 
tion rates), intrinsic noise levels in protein expression decrease. 
While the parameterisation of our dual reporter studies is such 
that these effects are small, the fact that less noise is expected 
at higher protein expression levels suggests a third mechanism 
by which high mitochondrial content may stabilise pluripotent 
cells. The contrasting results highlight the potential of exper- 
imental investigation of the effects of global transcription rate 
on the stability of multipotent states to inform of additional 
qualitative behaviors that models of lineage decision should be 
expected to exhibit. 



Discussion 

We have introduced a crude mathematical model for the 
effects of stochasticity in mitochondrial segregation and func- 
tionality on transcription rate in cells. Our model, while simple 
enough to allow some analytic treatment, reproduces a good 
number of experimentally observed features concerning the in- 
terplay of mitochondrial properties and transcription rate. We 
analyse our model and find that mitochondria provide extrin- 
sic noise contributions to transcription both through their un- 
even segregation at mitosis and through variability in their 
functionality. 

We note that, in addition to requiring variability in the 
amount of mitochondrial mass, an adequate fit to our data 
required us to consider variability in the function of mitochon- 
dria. This connects with the wealth of recent experimental 
and theoretical interest regarding the causes and control of 
heterogeneity of mitochondrial function [34l [35l |42j |44] and 
strengthens the case for the broad physiological relevance of 
functional variability. 

We incorporate our results for mitochondrial-sourced extrin- 
sic noise into existing models for mRNA and protein produc- 
tion, and show that mitochondrial noise can lead to significant 
variability between cells in a population. We also suggest that 
transcriptional variability resulting from mitochondrial noise 
may affect stem cell differentation, and illustrate this result 
with an analysis of two recent regulatory network-based mod- 
els for stem cell differentiation. We find that the quantitative 
effect of transcription rate variability on stem cell differenti- 
ation depends on the architecture of the regulatory network 
under consideration. 

Several recent studies have investigated the interplay be- 
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tween other possible sources of extrinsic noise in various or- 
ganisms. Before concluding, we will discuss connections to this 
body of literature. The recent study by Huh and Paulsson [29] 
found that variability in protein levels due to uneven inher- 
itance at mitosis might explain a body of experimental data 
that was previously assumed to result from noise in the protein 
production process. A mathematical study by Rausenberger 
et al. [21] also investigated the effects of inheritance stochas- 
ticity on cellular noise. Our work bears significant parallels 
to these ideas, in that we postulate uneven inheritance of mi- 
tochondria to be a substantial contributing factor to noise in 
all cellular processes that require ATP, including the mecha- 
nisms of protein production. Our philosophy also mirrors part 
of the work of Huh and Paulsson in that our model consid- 
ers a subset of cellular properties (in our case, mitochondrial 
partitioning and functionality, and cell volume) to provide all 
stochastic influences, with all other cellular properties evolving 
deterministically. 

The possible role of ATP as the proxy through which mito- 
chondrial variability affects other cellular processes ties in with 
an early prediction of Raser and O'Shea [M] who suggested 
that the dominance of extrinsic noise in expression variability 
across a wide range of proteins could result from fluctuations 
in a factor that affects expression for all genes. ATP, being re- 
quired for the processes of transcription and translation, meets 
this criterion. Shahrezaei et al. [45] illustrate the fact that ex- 
trinsic noise can influence intrinsic noise, through the former's 
effects on the rate constants involved in the latter. This in- 
fluence plays an important role in our model, where extrinsic 
variability of mitochondrial properties influences the synthe- 
sis rates of mRNA and protein through their dependence on 
[ylTP] . The ubiquity of ATP as an energy currency within the 
cell suggests that the rates of other intrinsic processes may be 
affected by the extrinsic variability we describe. 

The link between the process of transcription and noise in 
protein expression levels that we explore in the last section of 
this paper is related to the findings of Blake et al. [5? who 
found that protein expression noise depends on transcription 
efficiency. In our model, the modulation of transcription rate 
by noisy [ATP] has downstream effects on protein noise levels. 

Sigal et al. |26j , in a study of expression levels over a range of 
proteins, find cell cycle stage to be a significant contributor to 
extrinsic noise in protein abundance. Volfson et al. 19J used a 
mathematical framework to similarly identify population dy- 
namics, and upstream transcription factors, as key extrinsic 
contributors to cellular noise. Our model is compatible with 
these results, as cells at different cell cycle stages will have 
had different protein expression histories over their lifetimes. 
However, we anticipate that mitochondrial variability will also 
provide a significant contribution to protein expression noise, 
through modulation of upstream processes. 

An in-depth study by Newman et al. in yeast cells [TS] 
found a variety of protein-specific differences in expression 
noise according to transcription mode and protein function. 
Our model does not capture protein expression noise in this 
level of detail. The study of Newman et al. also characterised 
the contribution of intrinsic and extrinsic factors to total noise 
levels as a function of protein abundance. They found that 
while total expression noise did not scale with protein abun- 
dance, noise levels decreased with abundance when extrinsic 
factors were controlled for: suggesting that extrinsic factors 
were responsible for maintaining total noise levels as abun- 



dance increased. This suggestion that extrinsic noise increases 
in strength with protein abundance is captured in our protein 
level simulations. 

A study by Bar-Even et al., also in yeast cells [57], found 
intrinsic noise to be a substantial contributor to total noise, 
especially for proteins at intermediate abundance levels, with 
intrinsic contributions becoming less significant as expression 
levels increase (a similar result to Newman et al.). In this 
and several of the other studies above [13 [21], fluctuations 
in mRNA number were postulated to be the most impor- 
tant source of noise in protein expression levels. One mRNA 
molecule may produce many copies of a protein, and regulatory 
and chromatin-remodelling influences result in stochastic pro- 
duction of mRNA molecules, so random 'burstiness' in mRNA 
levels is a powerful source of noise in protein expression. 

Raj et al. pTfjj studied noise in mRNA expression in detail, 
and identify intrinsic effects as the dominant factors. Their 
study found that genes located in close proximity to each 
other displayed synchronised expression, while the expression 
of genes that were physically separate was unsynchronised, 
suggesting that local rather than global effects determine the 
expression levels of genes. While this study demonstrated that 
intrinsic effects significantly contribute to total noise in some 
cases, it was not explicitly shown that the magnitude of these 
effects outweighed extrinsic effects. Our results are compati- 
ble with this view that intrinsic noise plays an important role 
in gene expression, but we suggest that extrinsic noise due 
to energy variability may also be an important contributor to 
overall noise levels. 

We do not attempt to capture these mRNA processes ex- 
plicitly: rather, we take transcription rate to be a function of 
[ATP] as found in experiments [31] . However, we note the 
result that the measured functional form of this relationship 
changes in experiments in which chromatin was decondensed. 
This result suggests that the functional form of transcription 
rate with [ATP] allows us to capture some effects of the ATP- 
dependent chromatin remodelling process. 



Conclusions 

We find, through a phenomenological model constructed to 
reproduce recent data on mitochondrial and ATP variability, 
that stochastic inheritance of mitochondria at mitosis and vari- 
ability in mitochondrial function may be important sources 
of noise in transcription. By extension, these factors may 
contribute significantly to noise in protein expression further 
downstream. We have proposed experimental tests to refine 
our model and demonstrate its application in existing models 
for mRNA and protein production and stem cell differentation, 
and discussed how these findings integrate into the current un- 
derstanding of extrinsic noise in cellular biology. In particular, 
what our paper suggests is the need for multimodal single cell 
experiments through time (and through division) investigating 
coarse-grained measures of energy status, cellular volume, mi- 
tochondrial mass, and global rates of transcription and transla- 
tion. Cellular variability is of central physiological importance 
but we suggest that to understand this we must elucidate the 
relationships between certain core variables, including the rela- 
tionship between the machinery of expression and degradation 
and the energy status of the cell. 
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Methods 



Parameterisation of our model from experimental data. We 

used a subset of available experimental data (shown in Fig. [2l to choose 
numerical values for our key parameters. Some parameter values were 
fixed, in the sense that maximising the agreement between predictions 
from our model and a single experimental study allowed an optimal value 
to be chosen straightforwardly (si , c^ , and the ratio /3/a were parame- 
terised by data from Ref. 1311 (respectively the sigmoidal relationship be- 
tween [ATP] and transcription rate, the variance associated with volume 
partitioning, and the variance associated with mitochondrial partition- 
ing), v* by the maximal cell volume observed in data from Ref. |52| . and 
7 by comparing the mean properties of simulated cells to mean ATP lev- 
els reported in a population of HeLa cells in Ref. 1601 ). Other parameters 
influenced a range of predictions and values for these parameters were 
chosen by optimising the fit to experimental data across the set of results 
that they influenced (see 'Fitting Other Parameters' in Supplementary 
Information). 

Table 11] summarises the parameters and values employed in our model. 

Model Specifics. In the following we provide details of our model. 
The key equations specifying the key dynamics of cells in our model are 
given by Eqns. |1]4| The coupled ODEs in our model admit an analytic 
solution by simple integration; 



vit) = t,o-h^(exp(/3/t)-l), 
n(t) = noexp(/9/t). 



(5) 
(6) 



The reader will note that a variety of models for mitochondrial and vol- 
ume growth will yield similar forms. In addition, other functional forms 
for the level of ATP may be suggested. A selection of these alterna- 
tive forms are explored in 'Other Models' in Supplementary Information. 
Within this model, at mitosis, 



Vl_ 



ni 



M 



^liViV 



(7) 



(8) 



where Af{fJ., u) is a normal distribution with mean /i and variance a^, 
and V2 = V — vi,n2 = n — ni determine the daughter volumes fi,ii2 
and mitochondrial masses ni,n2. The variances of the mitochondrial 
distribution is chosen to represent a binomial distribution {B{n,p) ~ 
Minp, yjnp{l — p)), with p = ^ i for high n). The variance of the volume 
distribution is chosen to match experimental data on volume partitioning. 

Functionality evolves through changes at mitosis events, which follow 
an AR(1) process. A daughter cell's functionality /^ is determined from 
its parent's functionality /^: 

f''=faf'' + fc+N{0,af) (9) 

Upon mitosis, both daughter cells inherit the same /^ value, drawn 
from Eqn. |11| Once chosen, a cell's functionality remains constant 
throughout one cell cycle. We can model treatment with anti- and pro- 
oxidants by respectively increasing and decreasing /c, raising or lowering 
the mean functionality of mitochondria. With aj = 0.3, there is a very 
small but finite chance that cells will inherit zero (or lower) functionality. 
To avoid this unphysical case, we impose a cutoff of 0.01 on mitochondrial 
functionality. Values under this cutoff are resampled from Eqn. |11[ An- 
other model, where / varies continuously within a cell cycle, is considered 
in 'Other Models' in Supplementary Information. 

Virtual Mitochondria. In our model, the standard deviation in mi- 
tochondrial mass at mitosis is a function of the number of virtual mito- 
chondria in the cell, it„ = \/ti from binomial partitioning. We can vary 

the number of virtual mitochondria in the cell while keeping the total 
mitochondrial volume constant, by changing the volume assigned to an 
individual virtual mitochondrion. Let an individual virtual mitochondria 
have volume ^ , and let there be N virtual mitochondria in a cell. The 
total mitochondrial volume is Vn, and the mean inherited mitochondrial 
volume is half of this. The standard deviation in virtual mitochondrial 
number from a binomial distribution of the virtual mitochondria is 
so the standard deviation of mitochondrial volume is ^a. A standard 

deviation of a^ then corresponds to N = — ^ virtual mitochondria per 
cell. 

Cell Dynamics Simulations. To simulate a population of cells, we 
used a simple Euler method to solve the dynamic equations. A population 



of Af = 10 000 cells was simulated. When mitosis occurred, a random cell 
from the population was chosen for replacement by the new cell. To mea- 
sure distributions, this procedure was continued until the distributions 
stabilised. To measure sister-sister and between-generation correlations, 
a list of relationships between cells was maintained, with sister pairs only 
sampled when both sisters underwent an entire cell cycle without replace- 
ment. We note that this removal of randomly-chosen cells may potentially 
introduce artefacts into the results, as the constant probability of cell re- 
moval means that the probability of a cell surviving to a certain age is 
decreasing. We checked our results with a different simulation protocol; 
running the system and allowing exponential growth up to 10^ cells, with 
no removal. The resulting distributions were indistinguishable from the 
first protocol, showing that the removal rate is low enough so that the 
statistics are not affected. 

mRNA &L Protein Expression Simulations. For the simple sys- 
tem illustrating mRNA expression levels alone, cells were simulated us- 
ing the modification of Shahrezaei et al.'s 1451 to the Gillespie simu- 
lation method i63j in two scenarios; one involving only intrinsic noise 
((T,i = ay = af = 0) and one involving extrinsic noise due to mitochon- 
drial mass, functionality, and cell volume variability. In both cases, the 
initial copy number of mRNA molecules was set to zero, to illustrate 
differences in the dynamics of the system; transcription started at the 
birth of the cell, and there were no effects from mRNA inheritance. In 
the intrinsic noise experiments, a population of cells were simulated from 
identical initial conditions, with their n, v, f values set to the means of 
those variables obtained from simulations. The variability in mRNA ex- 
pression levels in these cells was therefore solely due to intrinsic noise. For 
the extrinsic noise simulations, mRNA expression was simulated in the 
heterogeneous population of cells that resulted from dynamic simulation. 
An ensemble of 2 500 cells was analysed for both cases and the mRNA 
content at each timestep recorded. 

To simulate the more complicated systems and investigate inheritance 
effects, we coupled the Shahrezaei et al. simulation protocol with the 
simple ODE solver so that simulation of a population of cells growing 
and producing mRNA and protein involved the following algorithm. 1) 
Use the ODE solver to calculate a cell's time of mitosis and the time 
series of volume and transcription rate throughout the cell lifetime. 2) 
Use Shahrezaei et al.'s method to compute the time behaviour of mRNA 
and protein levels given these time series for production rates. 3) Create 
daughter cells with noisy partitioning of volume, mitochondrial mass (bi- 
nomial) and functionality (AR(1)), and mRNA and protein copy numbers 
(binomial). 

After Raj et al. 1161 . we employ the following model values for 
birth (A) and death (5) rates of mRNA (m) and protein (n): (Am) = 
0.065-1, (A„> ^ 0.007 s-^,(m = 7 X 10-5s-i,^„ = 1.1 x lO-^s"!. In 
the case of birth rates, we used these mean values to scale the X{[ATP]) 
curves from das Neves et al. I31| so that the mean [ATP] level observed 
in a population gave the mean (A) values from Raj et al. (see 'mRNA & 
Protein Levels' in Supplementary Information). We also ran experiments 
with the degradation rates increased 100-fold; ^^n = 7 X 10"^ s~^,f„ = 
1.1 X 10""^ s~^, to explore the behaviour of the system at lower expression 
levels. 

Experiments. CMXRos labelling was done according to manufac- 
turer (invitrogen) instructions, cells were incubated with the probe (75 
nM) for 15 min. at 370 and washed with warm PBS. 

For the dual dye experiments, cells were loaded simultaneously with 
CMXRos and MitoTracker Green FM dye (Molecular Probes) for 20 min 
and after a brief PBS rinse fixed in PBS with 4% paraformaldehyde. 

Cell cycle length experiments were done by growing cells in the presence 
of media alone or containing either the anti-oxidant Dithiothreitol (250 
microM DTT) or the pro-oxidant Diamide (50 microM). Cell cycle length 
was measured as the time interval between two mitotic events in a single 
cell, analysed by live cell imaging using the Cell IQ platform and image 
analysis software (Chipman Tech.). 

In the fiow cytometry experiments, tripsinised Hela cells were washed 
and incubated with 20 nM MitoTracker Green FM dye (Molecular Probes) 
in medium at 37 C for 15 min, and then washed. Cells were analyzed using 
fiow cytometry (Dako CyAn ADP). 

Master equations. Let us consider the master equation for the tran- 
scription process, describing the probability of observing the system with 
m mRNAs at time t: 



dPm 

at 



: A(t)Pm_i + C(m + l)Pm+i - (A(t) + (m)P„ 



(10) 



where Pm is the probability of observing m mRNAs at a given time t, 
A(t) is transcription rate, and ^ is an mRNA degradation rate. 

Using a linear approximation for A(t) (see 'Parameterisation of A(t)' 
in Supplementary Information), we can solve this by using a generating 
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Parameter 


Description 


Value 


Motivation 


fa 


/ memory term 

Sets mean functionality 
(control) 


0.5 
0.5 


Fit parameter - chosen to give a mean 
functionality of 1 

Fit parameter - chosen to give a mean 
functionality of 1 


V* 

7 

Sl,2,3,4 


Volume for mitosis (scale) 

Proportionality between ^ 
and [ATP] 

V standard deviation at mi- 
tosis 

Set mean functionality 
(with anti-oxidant and 
pro-oxidant respectively) 
Fitting parameters for rela- 
tionship between [ATP] and 
A 


2 500 ^^m^ 
39 000 ^iM^m^ 
90 /im^ 
(0.69, 0.09) 

51.2, 44.7, 0.00288 ^M-\ 
-1.9 


Fixed for consistency with maximum vol- 
ume in Ref. |52| 

Fixed for consistency with mean ATP lev- 
els in Ref. [60] 

Fixed by volume segregation data in Ref. 
I31j 

Fixed by transcription rate noise levels in 
Ref. |31| 

Fixed by functional form of A in Ref. |3l| 


a 

/3 


V growth rate 

n growth rate 

/ standard deviation at mi- 
tosis 


0.92 hr--^ 

0.022 hr-^ 
0.34 


Chosen through optimisation - con- 
strained by mean cell cycle length in Ref. 
Ell- 
Fixed ratio with a through mitochondrial 
segregation data in Ref. J3lj 
Chosen through optimisation - con- 
strained through transcription rate noise 
and cell cvcle length variabilitv in Ref. J3l] 



TABLE I: Parameters and values employed in our 
Parameters' in Supplementary Information. 



model. For further information see 'Parameterisation of A(f)' and 'Fitting Other 



function approach (see 'mRNA & Protein Levels' in Supplementary In- 
formation). The mean, variance, and probability distribution of mRNA 
copy number at arbitrary time are given by: 






(11) 



(as + l)™o-le''i+"2 (ai + aiag + mo) 

(a3 + l)™»e«i+"^fai-faf + ^^ 
\ as + l 

xfl + 2ai + ^^^)-(a3 + l)"o+2 
V 0,3 + 1 J 

xe''i+'"2(ai-haia3-|-mo)2') (12) 

l_Fl(— m; mo — m -I- 1; — aias) (13) 



mo-m a2 
^3 



mo! 



m!(mo — m)! 
where iFi is the Kummer confluent hypergeometric function, a\ 



1. 



J [c + bt — ce ^' + |(e 5' — l)j , a2 = —moti; — ai and as = e^ 

Stem Cell Model. The progenitor cell differentation model of Huang 
et al. 1681 consists of the following equations for the evolution of protein 
expression levels xi (GATAl) and X2 (PU.l): 



dxi 
dt 

dx2 
dt 



+ x-i 



a2- 






- b2- 



kixi 

k2X2 



(14) 
(15) 



In this parameter set, a variables are self-activation rates, b vari- 
ables are cross-repression rates, k are decay rates, and 9 and n con- 
trol the functional form of these processes. Huang et al. show that 
altering these parameters changes the structure of the corresponding at- 
tractor landscape, so that the central undifferentiated attractor basin 
changes in size, affecting the predisposition of the system to differenti- 
ate. This landscape change gives a 'priming' of the system such that 
the effect of subsequent asymmetries may vary. In this study, we vary 
the landscape by symmetrically varying a and b, the parameters associ- 
ated with activation and repression, from the default parameterisation 
ai = bi = ki = 1, n = 4:, 9ai = d^i = 0.5 for i = 1,2. Specifically, we mod- 
ulated a and b with a multiplicative factor A; a — > Xa, b — > \b. We draw 



the connection between higher values of a and b and higher transcrip- 
tion and translation rates, as the rate of production of chemical species 
is increased by an increase in these parameters. 

The model of Chickarmane et al. |69) involves a relationshup between 
three dynamic variables: 



dt 

<m 

dt 

d[X] 
dt 



a-iA + a2[G] 
H-/3iA-f/32[G]-h/33[G][P] 
5iB + 52[P] 

■^^[G][P\^ 

-73[X] 



l + eiB + e2[P]- 

Cl[G] 
l + r,i[G] + mC 



-7i[G] 



-e4G][X] 



(16) 

-72[P] (17) 
(18) 



where [G], [P\, [X\ are the concentrations of GATAl, PU.l and a pos- 
tulated chemical species X respectively. A,B,C are external signalling 
factors: A and B promote GATAl and PU.l respectively, and C represses 
X. Other variables take default values oi = /3i = (5i = ei = /Ss = es = 1, 
"2 = /92 = ^2 = £2 = 0.25, £4 = 0.13, 71 = 72 = 73 = 0.01, 
Ci = ^1 = 0.01, 772 = 10. To vary transcription rates, we modulate 
a, 5 and C, terms with a multiplicative factor which we take to be pro- 
portional to transcription rate A. We used B = 0.5 with A = G = 0, 
allowing an external signal that promotes PU.l, to promote stability of 
the undifferentiated state. 
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FIG. 9: A with [ATP] (experiment and fit). 

Supplementary Information 
Parameterisation of A(t) 



Figs. 2g and 21 in das Neves et al. [31] give the response of transcription rate (A) in arbitrary units to varying concentrations 
of ATP, without and with artificial decondensation of chromatin respectively. Measuring [^TP] in iiM and working with the 
same arbitrary units employed in that study, we model this response with the expression: 



A = si + S2 taii-\s3[ATP] + S4), 



(19) 



with si ~ 51.2, S2 - 44.7, S3 ~ 2.88 x lO-^/zM"!, S4 ~ -1.9 (see Fig. |9|. 

das Neves et al. also produced a curve showing A with [ATP] in an experimental situation involving the decondensation 
of chromatin in the cell. This curve can be parameterised by the above equation with sf ~ 40.0, Sg — 54.7, S3 ~ 1.6 x 
10-^ fiM-\ si c^ -0.27. 

We choose this inverse tangent functional form to model the response of transcription rate to [ATP] for mathematical simplicity 
in subsequent sections, and note that modelling with other functional forms (for example. Hill functions) is also possible. 

In the parameterisation of our model, the time series of [ATP] in cells is rather linear and slowly varying. This behaviour 
emerges both due to the dynamics of mitochondrial density (which tends towards " with time) and the fact that the exponential 
growths involved are cell-cycle limited to only span a factor of around two before mitosis. It will be useful to employ a linear 
approximation to A(f), gained from an expansion in t about t': 



A(i) = si + S2 tan ^ S4 + S3 
~ c + bt, 
where the values of c and b are found after some algebra to be: 



7/7106' 



Pft 



vo+'-Ti^"'-^) 



(20) 
(21) 



b = 



S2f(3t'd 



(l + (s4 + rf)2) 

c — Si + S2 tan^"'"(s4 + d) 



and 



(22) 
(23) 



d = 



ano (e/3/*' - 1) + /Svo ' 



(24) 



Constants c and 6 then depend on cellular initial conditions and t' (which we take to be half the mean cell cycle length) but 
not on t. The sign of b is determined by the over- or under-population of mitochondria in the cell at mitosis: over-population 
will lead to high mitochondrial density and mean-reversion will act to decrease transcription rate with time (and vice versa). 
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Observable 


Experimental Value 


Simulated Value 


Error 


rjx Control 


0.4 


0.389 


0.028 


rjx Anti-oxidant 


0.2 


0.204 


0.019 


77a Pro-oxidant 


1 


0.943 


0.057 


r]\ Sister Cells 


0.08 


0.093 


0.166 


SD n+/n- 


0.23 


0.217 


0.055 


Cell cycle length 


30.5 


30.8 


0.009 


Mean [ATP] 


900 


925 


0.028 



TABLE II: Experimental observables and simulated fits obtained through optimisation of fitting parameters. Errors are the absolute 
difference between simulated and experimental values divided by experimental value. The score for a parameter set is the sum of errors 
for each observable. 
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FIG. 10: (A) Simulated relationship between a measure of total membrane potential in a cell (the product of mitochondrial mass n and 
functionality /) and the mitochondrial mass in the cell. (B) Membrane potential (measured with CMXRos) against mitochondrial content 
(measured with MitoGreen). 

Fitting Other Parameters 

The parameters in our model concerning volume cutoff and partitioning were chosen straightforwardly: v* for consistency 
with Ref. |52| and a^ to best fit experimental data from Ref. 31 . We then needed to find values for the remaining parameters 
influencing properties of our model cells, namely a and (3, the growth rates of cell volume and mitochondrial mass, 7, the constant 
of proportionality between — and [^TP], and parameters determining the statistics of mitochondrial functionality /a,/^,(T/. 
To fix values for these parameters, we first chose values for two parameters (fa and /^) to give a mean functionality of 1 in the 
control population, for simplicity. An optimisation procedure was then used to select values for af, a, /3 and 7. The procedure 
we used involved choosing arbitrary initial conditions for each parameter (though the order of magnitude of these initial values 
was determined by a crude preliminary investigation) and iterating, choosing a new value for one of the free parameters at each 
step, and retaining this value if the overall performance of the parameter set improved. The new values were chosen either 
(with probability 0.1) uniformly over the order of magnitude associated with that parameter or (with probability 0.9) uniformly 
from the interval of 10% deviations from the old value. To judge performance, a score for each parameter set was calculated 
based on the absolute deviation between experimental and simulated observables, averaged over the seven quantities shown in 
Table pi} Once these values were chosen for cells under 'normal' oxidative conditions (the default), f~'^'^ were set to match the 
experimentally observed transcription noise levels under the corresponding oxidative conditions. 



Mitochondrial Membrane Potential 



Fig. 



To] shows a comparison of the relationship between total membrane potential in a cell and mitochondrial mass in the 

Both our model and experimental data 
This result emerges straightforwardly 



cell, from new experiments (see Methods in Main Text) and simulation of our model. 

shows a linear correlation between total membrane potential and mitochondrial mass 

from our model due to the representation of total membrane potential (the product of n and /) and the qualitative agreement 

with experiment suggests that this modelling approach is suitable. 



Estimating Noise Contributions 



Here we find an expression for transcription rate in terms of the initial conditions and age of a cell, then estimate the variances 
associated with each of these contributory factors in order to assess the strength of each contribution to overall transcription 
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rate noise. 

As cells evolve deterministically within a cell cycle, the concentration of ATP in a model cell is a function of the cell's initial 
conditions and its age: 



[ATP] = 



7/n 



7/noe 



0ft 



vo + ^{ePf'-l)-' 



(25) 
(26) 



where the second line follows from simple integration of our dynamic equations. If we assume that the initial conditions of an 
individual cell are uncorrelated - which is the case for single divisions within our model, but population and memory effects may 
possibly lead to correlations - we can use standard rules for combining means and variances to calculate the population mean 
and variance of the ATP distribution from the total differential: 



fJ-ATP 



"ATP 



7M//^no 



gpt^ft^t 



M.o + ^(e^^^^'-l)' 



dATP 

dn„ 


2 

< + 


dATP 
df 


9 


dATP 

dvo 


< + 


dATP 
dt 



e2/t/3/32^2 



(nl {noa (e^*^ - 1 - ft/3) + v„/3{l + ftp)f a 



((e/*/5-l)noa + z;o/3)' 






(27) 
(28) 

(29) 
(30) 



As transcription rate A depends solely on [ATP] in our model, we can then calculate the mean and variance of the distribution 
of A: 



/^A 



^^ 



si 



d\ 



S2 tan 

2 



(S3^J'ATP + S4) 



2 
'^ATP 



dATP 

S2S3'^ATP 
1 + {iJ-ATPSz + 54)^ 



(31) 

(32) 

(33) 



From this, we can explore the contributions of mitochondrial segregation (c,^ ) and functional diversity {(j'^t) on transcription 



rate noise -qx 



^^. The overall expression can be written: 



■n\ = 



/"A 



S2S3CTATP 



(si + S2 tan ^{s^iiATP + S4))(l + (fJ-ATpSs + 54)^) 
S2S3jw?al^ + Wlal + W^ + W^aj 



(si + S2 tan ^{ssfiATP + Si)){^ + [tJ-ATpSz ■ 



^4))^ 



(34) 
(35) 

(36) 

(37) 



where the last line condenses the expression into the quadrature sum of noise levels in cellular parameters with weighting 
factors Wi — W^fif. 

Initial volume distribution. The mean and variance of initial volumes is simple to obtain, as division always occurs at a 
fixed volume and proceeds according to a fixed distribution. We straightforwardly have fi^ 



2 ' '^"o 




Initial mitochondrial mass distribution. We estimate the population variance of initial mitochondrial mass as the variance 
associated with the daughter of a cell with population average properties just before mitosis. Such a cell has v = v* 

The mean and variance associated with initial mitochondrial mass in the daughter of such a cell is fino = £ "T • 
These values give a simple estimate comparing to numerical results on the population statistics (see Fig. [llK). 

Mitochondrial functionality distribution. Similarly, we estimate the population distribution of mitochondrial function- 
alities by considering the daughter of a population-average cell. The steady-state distribution of the AR(1) process determining 

2 

/ is known to be normal with mean jzrr and variance y^W. It is not straightforwardly obvious that the distribution of / in 
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FIG. 11: A. Probability distribution of no under default model parameters, compared to the estimated distribution from considering a 
population-average cell. B. Probability distribution of / under default model parameters, compared to the distribution of the underlying 
AR(1) process. 
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FIG. 12: Probability distributions of cellular age t and cell cycle length f* under default model parameters. 

an unsynchronised population of cells, where the lifespan of a cell is a function of /, should also follow this form, but numerical 
results confirm that the / distribution does closely match it (see Fig. [Tl^). 

Distribution of cellular ages. This distribution is hard to estimate or approach analytically. We obtained the required 
values numerically from simulations, as shown in Fig. |12[ The mean and variance of this distribution was found to vary with 
cr„„ and af, the variance in mitochondrial mass and mitochondrial function. Over the ranges < cr„„ < 20 and < cr/ < 0.4, /it 
varied between 10 and 25 hours and at varied between 9 and 50 hours. However, these terms provide a negligible contribution 
to the transcription rate variability (see later), so we approximate the mean and variance of P(t) to be constant, with fit — 16 
and at — 14, the values corresponding to the default model parameterisation. 

We find, with our default parameter set, Wf = 8.5 x 10^,w„q = w^g — 4.2 x lO'^, and that Wt is zero within working 
precision. Typical noise levels in these quantities were measured from simulations as 77/ — 0.34, ?7„(j = 0.13, 77^^ = 0.07. This 
approximate analytic approach thus supports the hypothesis that the dominant contributions to transcription rate variability 
are from mitochondrial variability terms. 



Other Models 



Continuous f 



The behaviour of mitochondrial functionality / throughout and between cell cycles is the least well characterised property of 
our model. We use a simple picture in which / stays constant throughout a cell cycle and changes stochastically at mitosis, chosen 
for simplicity and due to the observed slowly varying behaviour of membrane potential in the cell cycle [21]. An alternative 
picture involves / being allowed to vary continuously and randomly throughout the cell cycle and no discrete jump occuring at 
mitosis, with daughter cells straightforwardly inheriting their parent's functionality. We found that this system was difficult to 
approach analytically, but performed numerical experiments to explore its behaviour. Instead of parameterising the behaviour 
of / with fa and /c, describing an autoregressive process, we now allow / to vary in a random walk with steps of size Sf per 
unit time. Sf then measures the degree of variability associated with /. To avoid unphysical situations, we place the constraint 
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FIG. 13; Time series with varying /. 
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FIG. 14: Behaviour of our model with continuously varying /. Comparison to experimental data of the model with / continuously 
varying within cell cycles and inherited continuously at mitosis. Plots are the same as Fig. [3] in the Main Text. 



0.1 < / < 1.9, forcing / to follow a constrained random walk. 

Figs. |13| and [l4| shows the behaviour of this model after parameterisation on the same data used to parameterise our default 
model. Most experimental results are captured to a similar extent than observed with our default model. These results suggest 
that the model's agreement with experimental data is robust with respect to the fine detail of the time evolution of / within cell 
cycles. This robustness suggests that it is the magnitude of extrinsic variability in mitochondrial functionality that gives rise to 
important physiological effects, rather than the specific detail of the time evolution and inheritance of functionality. However, as 
different dynamic regimes may adequately represent the behaviour of mitochondrial functionality, the need for further elucidation 
of mitochondrial functionality in model building is emphasised. 
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Label 


Equations 


Description 


A 


V — an 
h — j3n 


Default model type 


B 


V — av 
h — j3v 


Volume control on growth 


C 


V — an 

h =an{'y~^) 


Strong, explicit mean reversion of n towards a density 7. 


D 


V = Q« (l — :^) 
■h = Qn(7- ^) 


Explicit mean reversion of n towards 7, and volume 
growth towards a target v* 


E 


V — a 


Linear volume growth 


F 


i} — a 


Linear volume growth and weaker mean reversion on n 



TABLE III: Different possible models for time evolution of the cell exhibiting mean reversion. Each of these expressions admits an analytic 
solution for nif) and v(i). Other combinations (for example, weak mean reversion on n and exponential volume growth) are harder to 
solve and are often unstable for high or low initial njv. 

ATP and Alternatives 

It is thought that ATP levels in the cell are subject to homeostasis. For this reason we do not include a sink term for ATP, 
assuming that an ATP deficit will be immediately compensated for. A model that would capture ATP usage is given by: 



- /3' 



yn-{a' + f3')- 



(38) 
(39) 
(40) 



where a is the cellular ATP level. Here, the second term in Eqn. |40] corresponds to the rate of use of ATP in increasing 
the cell volume and mitochondrial mass. Note that if a = 772 our model equations are recovered. This model, with a suitable 
parameterisation, produces very similar results to our simpler model. 

We also note that the following system explicitly incorporates ATP homeostasis, as ATP is produced up to a certain level by 
mitochondria and used up in the increase of cellular volume: 



d = a'{a* — a)n — P'av, 
where a* controls the homeostatic level of ATP. This equation is solved by 



(41) 



aoe 



-t{a' n—fi' v) I 



a'n + 13' V ' 



(42) 



giving a hyperbolic form for ATP levels in terms of n and v. This representation could easily be extended to allow consideration 
of the ATP:ADP ratio ~ in which the rate of use of ATP gives the rate of production of ADP and vice versa. 

It is also possible that ROS plays a role in modulating transcription rate and cellular growth rates. In this case, an additional 
term could be introduced into the model, proportional to the number of mitochondria, and possibly inversely proportional to 
mitochondrial functionality, and used to modulate the dynamic equations and transcription rate. 

In our current model, we do not consider these complicating factors, due to our lack of experimental justification for them and 
also due to our desire to retain a simple, analytically tractable model. Future work, possibly motivated by the experiments we 
suggest in the main text, could take these more complicated factors into account. 

Other Deterministic Dynamic Forms 



Our model has been chosen phenomenologically to match experimental results concerning the distribution and behaviour of 
mitochondria. Other models are possible that display similar behaviour. Here we mention some potential alternative models. 
These alternatives both incorporate mean reversion on mitochondrial density and yield sensible key results. There is a difference 
in mean reversion rate between these models, illustrated by dynamic plots in Fig. [15] 

Our model was chosen for its analytic tractability and its ability to reproduce experimental phenomena: these more complex 
models rapidly become analytically intractable, which hinders a deeper understanding of their behaviour. 
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FIG. 15: The functional form of mitochondrial density p with time t (in arbitrary units) in other possible models for the time evolution 
of cellular properties. For these illustrations, the parameter set a — 0.1, /3 = 0.01, fo = 1000 is used. The different lines correspond to 
trajectories of p resulting from varying initial no between and 400 (a much wider range than in the default model, as different choices of 



functional form may lead to different absolute values for n and v). The label of the plot gives the corresponding function in Table III 



We note in addition that the [ATP] term that appears in our dynamic equations may be replaced by a transcription or 
translation rate term, explicitly including the sigmoidal (or hyperbolic) response of these cellular growth rates in the dynamic 
equations. In our bare model, we use the simpler [ATP] term both to maintain analytic tractability and because the appropriate 
form for a response curve is difficult to estimate (for example, cellular growth may be a result of a combination of transcription 
and translation, and may include other additional terms). This approach essentially involves using a linear approximation for 
growth rate: linear approximations for transcription rate are employed elsewhere in this study and yield very similar results to 
the numerical behaviour of the corresponding nonlinear case. 

It is possible to construct models for mitochondrial stochasticity that do not involve mean-reverting behaviour of mitochondrial 
density. As a simple example, if cellular dynamics are such that mitochondrial mass and cellular volume evolve independently: 



V 

n 



afv; 
Pfn, 



(43) 

(44) 



with mitosis involving a binomial segregation of volume and mitochondrial mass as before, the population variance in mitochon- 
drial density will generally increase in an unbounded manner, as the probability of observing higher and higher mitochondrial 
densities increases with time. However, if cutoffs are placed on mitochondrial density, so that cells are removed from the 
population if their density does not obey 



n 
<- <P+, 

V 



(45) 



the population distribution of mitochondrial density can achieve stationarity. This approach results in a much weaker cor- 
relation between mitochondrial mass n and cellular volume v. This model can be parameterised (an example parameter set is 
p_ = 10, p+ = 300, a = /3 = 0.1,7 = 3.6, other parameters the same as in the default model) to yield similar results to the 
default model for transcriptional noise. 

This rather constructed model is considered due to the discrepancy between flow cytometry data in das Neves et al. j31) . 
showing an absence of correlation between mitochondrial volume and cellular volume. Obtaining analytic results from this model 
is difficult, both in terms of approximation for the moments of probability distributions and the time evolution of mRNA. 



Potential Experiments for Refinement 



As mentioned in the Main Text, our model was constructed from a phenomenological philosophy, with the intention of using 
experimental results to construct a plausible coarse-grained explanation for the influence of mitochondrial variability on extrinsic 
noise in general and transcription rate in particular. Our goal was to introduce a simplified but consistent mathematical summary 
of the data and to use this to motivate further experiments. We suggest a set of experiments in Table 11 that would support or 
contribute to further development of this model. 
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Observable 



Mitochondrial functionality (measured via 
membrane potential) before and after mi- 
tosis. 



Protein expression levels as a function of 
[ATP]. 

Time series of mitochondrial mass and den- 
sity through the cell cycle. 



Behaviour of [ATP] with time. 

Determinant factors of [ATP]. 

Relative contribution of mitochondrial 
mass and functional variability to tran- 
scriptional variability (measured via bro- 
mouridine uptake). 
Noise in protein expression levels. 

ROS levels (measured via, for example, Mi- 
toSox). 



Prediction 



Parent membrane potential is weakly retained but 
stochastically altered to give daughter membrane poten- 
tial. Note that our current model treats this process 
crudely and as such many other dynamics may be possi- 
ble. 

Expression levels are higher in cells with higher [ATP]. 

Exponential cell growth, with mitochondrial density 
tending towards an average value. The time series of 
these dynamics could be used to distinguish between dif- 
ferent models (see 'Other Models'). 

Slowly varying over the cell cycle, mean-reverting to- 
wards an average value. 

Mitochondrial mass and membrane potential, and cell 
volume. Levels of ROS may also be of importance. 
Strong dependence of transcriptional noise on both mi- 
tochondrial inheritance and functionality variability. 



Dependent on mitochondrial variability, and lower for 
cells with mitochondrial low mass and functionality. 
Possibly correlated inversely with a measure of mitochon- 
drial functionality, and possibly higher with low mito- 
chondrial mass and functionality, as weaker/sparser mi- 
tochondria struggle to match energy demands. 



TABLE IV: Potential experiments that may clarify aspects of our model, roughly ranked in order of importance for supporting or suggesting 
area of refinement for our model. 



mRNA & Protein Levels 



The master equation for the probability distribution of mRNA levels can be written down as: 



dPm 
dt 



A(t)P,„_i + e(m + l)P,n+i - (A(t) + ^m)P„„ 



(46) 



where A(t) is the time-dependent birth rate of mRNA molecules, ^ is the rate of removal of mRNA, and -P,„(t) is the probability 
of observing the system with m mRNA molecules at time t. We choose 



A(i) = (c + fet), 

the approximation from 'Parameterisation of A(t)', to model transcription rate with time. 
Using the generating function G{z) — J2m ^™^m gives: 



(47) 



- = A(t)(.-l)G-e(.-l)^. 



(48) 



We can solve Eqn. |48] with the method of characteristics. This process allows us to convert the PDE into a set of ODEs 



along a characteristic curve of the function. We wish to recast Eqn. 48 into the following form, where a characteristic curve is 
parameterised by the new variable s: 



Using the chain rule, we can write: 



ds 



G{z{s),t{s))=F{G,z{s),tis)). 



(49) 



ds 



G{zis),t{s)) 



dG dz dG dt 

dz ds dt ds 

1 dG BG 



z — 1 dt dz 



^—^{c + bt)G, 



(50) 
(51) 
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where the last hne is just a rearrangement of the original PDE. By comparing coefficients, we obtain the characteristic equations 
for the system: 



dz 

ds 
dt 

ds 

dG 

ds 



£, -^ z ^ ^s + Zq 
1 1 



Z - f ^S + Zq — f 



-> t 



ln($s + zo - 1) + io, 



e 



(c + bt)G =\c + b\- ln(es + zo - f ) + to 1 ) G. 



(52) 
(53) 
(54) 



These equations may be interpreted as describing the evolution of the arguments of G, and the subsequent evolution of the 
value of G, as we progress along a curve parameterised by s. As absolute values of s are not important, affecting only the 
parameterisation of progress along a characteristic curve, we set zq — Q without loss of generality. We then have s = f and 
io = t — I ln(z — 1). The final ODE can then be solved by separation of variables: 



G 



-dC = 



b ( - ln(^s + zo - 1) + C2 ) ]ds 



s(c + Wo) - -J 

z , , , bz 
^(c + 5t)--- 



sb , , ^ , b , , 

- HCs - 1) - ^ ln(Cs 

^ln(z-l), 



(55) 
(56) 
(57) 



We then have 



G-Gexp(|(c + 6t)-pK^-l)^' 



(58) 



where the arbitrary function C = C{t — j ln(z — 1)) = G(to), as this quantity is independent of s. 

If the initial copy number of mRNA molecules is mo, we have the initial condition G{z, 0) = J2m ^™Prn{0) — J2m ^"^^n 
"^° . Noting that at i = 0, 6"*"^ = z — 1, we find that if we employ the choice 



C(to) - exp (^{e-'oi + 1) (^^ _ ^^ j (e-*o« + i)m„ e^p (^-^to(^ 



(59) 



we recover the required initial condition: 



bz 



-to? + lyno g^p ,^^ 



-b 



G(z,to) = exp(^-(c + 6i)--j(z-l)5^exp(^(e-*°« + l)^^2 ^JJ- —' -- y^2 

G(z,to = 0) = exp(^|-p^(z-l)i^cxp(^z(^A-^)^z™°(z-l)^ 



(60) 

(61) 
(62) 



The general solution is then given by Eqn. 60 which, using tg — t — j ln(z — 1), can be written: 



G(z,i) = e"i"+'^H^ + a3) 



mo 



(63) 



with 



ai 

02 

as 



bt 



<t 



+ ^(e-^*-l)) 



— moi^ — oi 



(64) 

(65) 
(66) 



We can recover the mean and variance of the distribution Pm{t) by using ^„ 



WL=i anda2„ 






dG l 

dz \z=l 



dG l ^2. 
dz \z=l) ■ 



(a.3 + l)'"«-ie''i+°= (ai + a^a^ + mo) 
\ aa + l 



03 + 1 
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(67) 
(68) 



We can also recover the full probability distribution of observing m mRNAs using Pm{i) — ^jb^z^l -o' Using Leibniz's rule 
allows us to write: 



Pm{t) 



1 a^G 



m\ dz 
1 



E 



2=0 

TO! •' — ' \ I J dz' 
=0 



y -T, — ^«ie" 

-^-^zKto-i)! 






c^z"^ 



2=0 



Too! 



12+02 "" • ('^ I „ \ma-m+i 

(ttio — m + i)\ 



2=0 



E 

i=0 



molale^ftg 
(too — to + i)!i!(TO — i)! 



This can be alternatively written as 



Pm{t) 



Too! 



?7l! (777,9 



iFi(— to;too — 777 + 1; —0103) 



where iFi is the Kummcr confluent hypergeometric function. 
The corresponding master equation including proteins is: 



(69) 
(70) 

(71) 
(72) 



(73) 



^^ rn n 

dt 



= m\n{t)Pjnn-l + Cn{n + l)P„„+i - £,nnPmn " TOA„(t)P„i„ 



(74) 



for which we have not obtained a general analytic solution. Instead, we simulate this system with a stochastic simulation algo- 
rithm. To simulate these systems, we use a parameter set employed by Raj et al. |16) : (Am) = 0.06s~^, (A„) = 0.007s~^,^„i = 
7 X 10^^ s^^, ^„ = 1.1 X 10~^ s~^. The death rates are simply inserted into the simulation, and the mean birth rates are used to 
tune the time- varying birth rate for the simulation. For example, the normalisation of Xmit) was chosen so that the population 
average value of [ATP] corresponded to (Am)- 

As the proces s of translation is believed to be [ATPj-dependent but not dependent on chromatin remodelling, the decondensed 
version of Eqn. 19 (using the hyperbolic sf coefficients rather than the sigmoidal Si coefficients) was used for A„. Again, the 



mean birth rate was used to normalise this function as above. Overall, we then have constant rates for mRNA death (faster) and 
protein death (slower) and time-varying birth rates for mRNA (including chromatin remodelling) and proteins (not including 
chromatin remodelling). 

As mentioned in the Main Text, we found this parameterisation to yield large copy numbers of proteins, which led to very 
low values for intrinsic noise. We explored this effect by increasing the degradation rates from Raj et al.'s default values, to 
^m = 7 X 10^'^ s^^iCn — 1-1 X 10"'^. In addition, we investigated the case where variability in [ATP] only affected transcription 
rate, with translation rate taking the constant value A„ = 0.007 s""'^. As shown in the Main Text, those simulations with higher 
degradation rates showed a more significant intrinsic noise contribution, and the trends in these dual reporter simulations changed 
little with a constant translation rate. This consistency means that though our crude model of the dependence of translation 
rate on [ATP] could be challenged, even if translation is ATP dependent we still see pronounced effects on protein levels through 
ATP-dependent variability in transcription rate. 
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